Wisconsin Prognosis

Libraries

library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
#pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)

The data

dataBreast <- read.csv("~/GitHub/RISKPLOTS/DATA/wpbc.data", header=FALSE)
table(dataBreast$V2)
## 
##   N   R 
## 151  47
rownames(dataBreast) <- dataBreast$V1
dataBreast$V1 <- NULL
dataBreast$status <- 1*(dataBreast$V2=="R")
dataBreast$V2 <- NULL
dataBreast$time <- dataBreast$V3
dataBreast$V3 <- NULL
dataBreast <- sapply(dataBreast,as.numeric)
## Warning in lapply(X = X, FUN = FUN, ...): NAs introduced by coercion
dataBreast <- as.data.frame(dataBreast[complete.cases(dataBreast),])
table(dataBreast$status)
## 
##   0   1 
## 148  46

Modeling

ml <- BSWiMS.model(Surv(time,status)~1,data=dataBreast)

[++++]

sm <- summary(ml)
pander::pander(sm$coefficients)
Table continues below
  Estimate lower HR upper u.Accuracy r.Accuracy
V26 8.07e-03 1 1.01 1.01 0.593 0.237
V27 4.13e-04 1 1.00 1.00 0.608 0.237
V24 7.71e-03 1 1.01 1.01 0.598 0.634
V35 8.65e-06 1 1.00 1.00 0.727 0.237
V34 9.13e-03 1 1.01 1.02 0.634 0.598
Table continues below
  full.Accuracy u.AUC r.AUC full.AUC IDI NRI z.IDI
V26 0.593 0.598 0.500 0.598 0.0626 0.393 2.77
V27 0.608 0.608 0.500 0.608 0.0563 0.434 2.76
V24 0.603 0.609 0.618 0.613 0.0532 0.323 2.62
V35 0.727 0.641 0.500 0.641 0.0289 0.565 2.28
V34 0.603 0.618 0.609 0.613 0.0233 0.411 2.13
  z.NRI Delta.AUC Frequency
V26 2.38 0.09827 1
V27 2.63 0.10840 1
V24 1.94 -0.00529 1
V35 3.50 0.14116 1
V34 2.47 0.00338 1

Cox Model Performance

Here we evaluate the model using the RRPlot() function.

The evaluation of the raw Cox model with RRPlot()

Here we will use the predicted event probability assuming a baseline hazard for events withing 5 years

index <- predict(ml,dataBreast)
timeinterval <- 2*mean(subset(dataBreast,status==1)$time)

h0 <- sum(dataBreast$status & dataBreast$time <= timeinterval)
h0 <- h0/sum((dataBreast$time > timeinterval) | (dataBreast$status==1))
pander::pander(t(c(h0=h0,timeinterval=timeinterval)),caption="Initial Parameters")
Initial Parameters
h0 timeinterval
0.323 51.1 
rdata <- cbind(dataBreast$status,ppoisGzero(index,h0))
rownames(rdata) <- rownames(dataBreast)

rrAnalysisTrain <- RRPlot(rdata,atProb=c(0.90),
                     timetoEvent=dataBreast$time,
                     title="Raw Train: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

As we can see the Observed probability as well as the Time vs. Events are not calibrated.

Uncalibrated Performance Report

pander::pander(t(rrAnalysisTrain$OERatio),caption="O/E Ratio")
O/E Ratio
est lower upper
0.843 0.617 1.12 
pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Ratio")
O/E Ratio
mean 50% 2.5% 97.5%
1 1 0.955 1.05 
pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Ratio")
O/Acum Ratio
mean 50% 2.5% 97.5%
0.776 0.775 0.768 0.783 
pander::pander(rrAnalysisTrain$c.index$cstatCI,caption="C. Index")
mean.C Index median lower upper
0.677 0.677 0.595 0.755 
pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.635 0.543 0.727 
pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.261 0.143 0.411 
pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.899 0.838 0.942 
pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90%
0.419 
pander::pander(t(rrAnalysisTrain$RR_atP),caption="Risk Ratio")
Risk Ratio
est lower upper
1.88 1.1 3.2 
pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")
Logrank test Chisq = 11.670372 on 1 degrees of freedom, p = 0.000635
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 167 34 41.11 1.23 11.7
class=1 27 12 4.89 10.33 11.7

Cox Calibration

op <- par(no.readonly = TRUE)


calprob <- CoxRiskCalibration(ml,dataBreast,"status","time")

pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Cox Calibration Parameters")
h0 Gain DeltaTime
0.248 0.767 41

The RRplot() of the calibrated model

h0 <- calprob$h0
timeinterval <- calprob$timeInterval;

rdata <- cbind(dataBreast$status,calprob$prob)


rrAnalysisTrain <- RRPlot(rdata,atProb=c(0.90),
                     timetoEvent=dataBreast$time,
                     title="Calibrated Train: Breast",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

Calibrated Train Performance

pander::pander(t(rrAnalysisTrain$OERatio),caption="O/E Ratio")
O/E Ratio
est lower upper
0.877 0.642 1.17 
pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Ratio")
O/E Ratio
mean 50% 2.5% 97.5%
1.05 1.05 1 1.1 
pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Ratio")
O/Acum Ratio
mean 50% 2.5% 97.5%
0.943 0.943 0.933 0.953 
pander::pander(t(rrAnalysisTrain$c.index$cstatCI),caption="C. Index")
C. Index
mean.C Index median lower upper
0.677 0.678 0.595 0.756 
pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.635 0.543 0.727 
pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.261 0.143 0.411 
pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.899 0.838 0.942 
pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90%
0.341 
pander::pander(t(rrAnalysisTrain$RR_atP),caption="Risk Ratio")
Risk Ratio
est lower upper
1.88 1.1 3.2 
pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")
Logrank test Chisq = 11.670372 on 1 degrees of freedom, p = 0.000635
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 167 34 41.11 1.23 11.7
class=1 27 12 4.89 10.33 11.7

Cross-Validation

Here we use the estimated h0 and timeinterval from the full set

rcv <- randomCV(theData=dataBreast,
                theOutcome = Surv(time,status)~1,
                fittingFunction=BSWiMS.model, 
                trainFraction = 0.9,
                repetitions=100,
                classSamplingType = "Pro"
         )

.[+++++].[+++].[++++].[++++].[-+].[–].[–].[++++++].[+++++].[++++++]10 Tested: 118 Avg. Selected: 3.4 Min Tests: 1 Max Tests: 5 Mean Tests: 1.694915 . MAD: 0.452554

.[+++].[+++].[++++++].[++++].[+++++].[+++++].[++++++].[+++++++++].[++].[++++]20 Tested: 167 Avg. Selected: 4.35 Min Tests: 1 Max Tests: 6 Mean Tests: 2.39521 . MAD: 0.4697117

.[+++].[+++++].[+++++].[+++++].[+++++].[+++++++].[+++++++].[++++++].[+++].[+++++++]30 Tested: 184 Avg. Selected: 4.866667 Min Tests: 1 Max Tests: 8 Mean Tests: 3.26087 . MAD: 0.4737094

.[++++++].[++++].[+++].[++++++].[+++++].[++++].[++++].[–].[++++].[++++++]40 Tested: 192 Avg. Selected: 4.925 Min Tests: 1 Max Tests: 10 Mean Tests: 4.166667 . MAD: 0.4739559

.[+++++].[+++++++].[++++++].[++++].[+++++].[++++++].[++++++].[++++++].[++++].[++++]50 Tested: 193 Avg. Selected: 5.12 Min Tests: 1 Max Tests: 11 Mean Tests: 5.181347 . MAD: 0.4751408

.[+++].[++++++++].[++++++].[++++].[+].[++].[++++++++].[++++++++].[++++++].[+++]60 Tested: 194 Avg. Selected: 5.183333 Min Tests: 1 Max Tests: 13 Mean Tests: 6.185567 . MAD: 0.4779008

.[-+++].[++++].[++++].[+++].[++++++].[+++++++++].[++++++].[+++++].[++++].[++++]70 Tested: 194 Avg. Selected: 5.185714 Min Tests: 1 Max Tests: 15 Mean Tests: 7.216495 . MAD: 0.479849

.[+++].[++++].[–].[+++].[+++++].[++++++].[+++].[+++++++].[+++++++].[++]80 Tested: 194 Avg. Selected: 5.0875 Min Tests: 1 Max Tests: 17 Mean Tests: 8.247423 . MAD: 0.4792785

.[++++].[+++++].[+++++].[++++].[++++].[+++].[+++++].[++++].[+].[+++++++]90 Tested: 194 Avg. Selected: 5.033333 Min Tests: 2 Max Tests: 19 Mean Tests: 9.278351 . MAD: 0.4819606

.[++++].[++++++++].[++++].[+++].[+].[++++].[+].[++++++].[+].[+++++++]100 Tested: 194 Avg. Selected: 4.95 Min Tests: 2 Max Tests: 20 Mean Tests: 10.30928 . MAD: 0.4823863

stp <- rcv$survTestPredictions
stp <- stp[!is.na(stp[,4]),]

bbx <- boxplot(unlist(stp[,1])~rownames(stp),plot=FALSE)
times <- bbx$stats[3,]
status <- boxplot(unlist(stp[,2])~rownames(stp),plot=FALSE)$stats[3,]
prob <- ppoisGzero(boxplot(unlist(stp[,4])~rownames(stp),plot=FALSE)$stats[3,],h0)

rdatacv <- cbind(status,prob)
rownames(rdatacv) <- bbx$names
names(times) <- bbx$names

rrAnalysisTest <- RRPlot(rdatacv,atThr = rrAnalysisTrain$thr_atP,
                     timetoEvent=times,
                     title="Test: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

Cross-Validation Test Performance

pander::pander(t(rrAnalysisTest$OERatio),caption="O/E Ratio")
O/E Ratio
est lower upper
0.864 0.633 1.15 
pander::pander(t(rrAnalysisTest$OE95ci),caption="O/E Ratio")
O/E Ratio
mean 50% 2.5% 97.5%
1.04 1.04 0.99 1.09 
pander::pander(t(rrAnalysisTest$OAcum95ci),caption="O/Acum Ratio")
O/Acum Ratio
mean 50% 2.5% 97.5%
0.904 0.904 0.891 0.916 
pander::pander(rrAnalysisTest$c.index$cstatCI,caption="C. Index")
mean.C Index median lower upper
0.655 0.656 0.57 0.734 
pander::pander(t(rrAnalysisTest$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.607 0.512 0.702 
pander::pander((rrAnalysisTest$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.239 0.126 0.388 
pander::pander((rrAnalysisTest$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.865 0.799 0.915 
pander::pander(t(rrAnalysisTest$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90%
0.341 
pander::pander(t(rrAnalysisTest$RR_atP),caption="Risk Ratio")
Risk Ratio
est lower upper
1.55 0.884 2.72 
pander::pander(rrAnalysisTest$surdif,caption="Logrank test")
Logrank test Chisq = 5.216087 on 1 degrees of freedom, p = 0.022379
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 163 35 40.14 0.657 5.22
class=1 31 11 5.86 4.498 5.22

Calibrating the test results

rdatacv <- cbind(status,prob,times)
calprob <- CalibrationProbPoissonRisk(rdatacv)

pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Cox Calibration Parameters")
h0 Gain DeltaTime
0.323 1 41.3 
timeinterval <- calprob$timeInterval;

rdata <- cbind(status,calprob$prob)


rrAnalysisTest <- RRPlot(rdata,atProb=c(0.90),
                     timetoEvent=times,
                     title="Calibrated Test: Breast",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

### Calibrated Test Performance

pander::pander(t(rrAnalysisTest$OERatio),caption="O/E Ratio")
O/E Ratio
est lower upper
0.871 0.638 1.16 
pander::pander(t(rrAnalysisTest$OE95ci),caption="O/E Ratio")
O/E Ratio
mean 50% 2.5% 97.5%
1.05 1.05 1 1.11 
pander::pander(t(rrAnalysisTest$OAcum95ci),caption="O/Acum Ratio")
O/Acum Ratio
mean 50% 2.5% 97.5%
0.904 0.903 0.891 0.916 
pander::pander(rrAnalysisTest$c.index$cstatCI,caption="C. Index")
mean.C Index median lower upper
0.655 0.653 0.566 0.737 
pander::pander(t(rrAnalysisTest$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.607 0.512 0.702 
pander::pander((rrAnalysisTest$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.196 0.0936 0.339 
pander::pander((rrAnalysisTest$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.899 0.838 0.942 
pander::pander(t(rrAnalysisTest$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90%
0.369 
pander::pander(t(rrAnalysisTest$RR_atP),caption="Risk Ratio")
Risk Ratio
est lower upper
1.68 0.931 3.04 
pander::pander(rrAnalysisTest$surdif,caption="Logrank test")
Logrank test Chisq = 4.701728 on 1 degrees of freedom, p = 0.030132
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 170 37 41.4 0.467 4.7
class=1 24 9 4.6 4.196 4.7